In [1]:
%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import pi, var, S
from sympy.utilities.codegen import codegen


Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

The following code is taken from the article:

Čertík, O., Winkler, P. (2013). Computation of screened two-electron matrix elements. International Journal of Quantum Chemistry. doi:10.1002/qua.24431


In [2]:
from sympy import var, legendre, Integral, \
    exp, latex
var("l R alpha t")
f = (2*l+1) / (2*t) * Integral(legendre(l, \
    (1-R**2+t**2) / (2*t)) * \
    exp(-alpha*R), (R, 1-t, 1+t))
for _l in range(3):
    expr = f.subs(l, _l).doit().simplify( ) \
        / exp(-alpha)
    expr = expr.series(alpha, 0, 14)
    print " s_%d =" % _l, latex(expr)


 s_0 = 1 + \frac{1}{6} \alpha^{2} t^{2} + \frac{1}{120} \alpha^{4} t^{4} + \frac{1}{5040} \alpha^{6} t^{6} + \frac{1}{362880} \alpha^{8} t^{8} + \frac{1}{39916800} \alpha^{10} t^{10} + \frac{1}{6227020800} \alpha^{12} t^{12} + \mathcal{O}\left(\alpha^{14}\right)
 s_1 = t + \alpha t + \frac{1}{10} \alpha^{2} t^{3} + \frac{1}{10} \alpha^{3} t^{3} + \frac{1}{280} \alpha^{4} t^{5} + \frac{1}{280} \alpha^{5} t^{5} + \frac{1}{15120} \alpha^{6} t^{7} + \frac{1}{15120} \alpha^{7} t^{7} + \frac{1}{1330560} \alpha^{8} t^{9} + \frac{1}{1330560} \alpha^{9} t^{9} + \frac{1}{172972800} \alpha^{10} t^{11} + \frac{1}{172972800} \alpha^{11} t^{11} + \frac{1}{31135104000} \alpha^{12} t^{13} + \frac{1}{31135104000} \alpha^{13} t^{13} + \mathcal{O}\left(\alpha^{14}\right)
 s_2 = t^{2} + \alpha t^{2} + \alpha^{2} \left(\frac{1}{14} t^{4} + \frac{1}{3} t^{2}\right) + \frac{1}{14} \alpha^{3} t^{4} + \alpha^{4} \left(\frac{1}{504} t^{6} + \frac{1}{42} t^{4}\right) + \frac{1}{504} \alpha^{5} t^{6} + \alpha^{6} \left(\frac{1}{33264} t^{8} + \frac{1}{1512} t^{6}\right) + \frac{1}{33264} \alpha^{7} t^{8} + \alpha^{8} \left(\frac{1}{3459456} t^{10} + \frac{1}{99792} t^{8}\right) + \frac{1}{3459456} \alpha^{9} t^{10} + \alpha^{10} \left(\frac{1}{518918400} t^{12} + \frac{1}{10378368} t^{10}\right) + \frac{1}{518918400} \alpha^{11} t^{12} + \alpha^{12} \left(\frac{1}{105859353600} t^{14} + \frac{1}{1556755200} t^{12}\right) + \frac{1}{105859353600} \alpha^{13} t^{14} + \mathcal{O}\left(\alpha^{14}\right)

In [3]:
f


Out[3]:
$$\frac{\left(2 l + 1\right) \int_{- t + 1}^{t + 1} e^{- R \alpha} P_{l}\left(\frac{- R^{2} + t^{2} + 1}{2 t}\right)\, dR}{2 t}$$

In [4]:
expr


Out[4]:
$$t^{2} + \alpha t^{2} + \alpha^{2} \left(\frac{1}{14} t^{4} + \frac{1}{3} t^{2}\right) + \frac{1}{14} \alpha^{3} t^{4} + \alpha^{4} \left(\frac{1}{504} t^{6} + \frac{1}{42} t^{4}\right) + \frac{1}{504} \alpha^{5} t^{6} + \alpha^{6} \left(\frac{1}{33264} t^{8} + \frac{1}{1512} t^{6}\right) + \frac{1}{33264} \alpha^{7} t^{8} + \alpha^{8} \left(\frac{1}{3459456} t^{10} + \frac{1}{99792} t^{8}\right) + \frac{1}{3459456} \alpha^{9} t^{10} + \alpha^{10} \left(\frac{1}{518918400} t^{12} + \frac{1}{10378368} t^{10}\right) + \frac{1}{518918400} \alpha^{11} t^{12} + \alpha^{12} \left(\frac{1}{105859353600} t^{14} + \frac{1}{1556755200} t^{12}\right) + \frac{1}{105859353600} \alpha^{13} t^{14} + \mathcal{O}\left(\alpha^{14}\right)$$

In [5]:
[(fname, fcode), (hfile, hcode)] = codegen(("f", expr), "F95", "test", header=False, empty=False)

In [6]:
fcode


Out[6]:
                       REAL*8 function f(alpha, t)                        
                              implicit none                               
                       REAL*8, intent(in) :: alpha                        
                         REAL*8, intent(in) :: t                          
                          REAL*8 :: O(alpha**14)                          
f = t**2 + alpha*t**2 + alpha**2*(t**4/14 + t**2/3) + alpha**3*t**4/14 + &
      alpha**4*(t**6/504 + t**4/42) + alpha**5*t**6/504 + alpha**6*(t** & 
        8/33264 + t**6/1512) + alpha**7*t**8/33264 + alpha**8*(t**10/ &   
      3459456 + t**8/99792) + alpha**9*t**10/3459456 + alpha**10*(t**12 & 
       /518918400 + t**10/10378368) + alpha**11*t**12/518918400 + alpha & 
       **12*(t**14/105859353600 + t**12/1556755200) + alpha**13*t**14/ &  
                          105859353600 + O(alpha**14)                     
                               end function                               

In [6]: